#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>

#include "Visualizator/visualizator.03.h"
#include "Vector/vector.03.h"
#include "Color/color.01.h"
#include "Timer/timer.01.h"
#include <vector>
#include <deque> 

using namespace std;

const double kB = 1.3806488e-23;

const int gsizex = 100;
const int gsizey = 100;
const int gsizez = 100;

double gstepx = 3;
double &gstepy = gstepx;
double &gstepz = gstepx;


double e = 127*12;	
double drij2max = SQR(3.0); 
double dt = 1e-3;

typedef double real;

struct P {
	Vector<double> R;
	Vector<double> V;
	Vector<real> F;
	int type;
};


vector<P*> g[gsizex][gsizey][gsizez];

class MyVisualizator: public Visualizator {
	public:

		struct PW {
		
			double e;
			double s; 
			double rc;
		
		};

		struct PT {		

			double m;
			Color color;	
		
		};
		
		
		char **argv;
		
		PW pw[10][10];
		
		vector< PT > pt;
		
		vector<P> p;
		

		
		bool coolingF;
		
		double rnd() {
			return (double)rand()/RAND_MAX;
		};
		
		
		void packToHCP(vector<P> &p, PW &pw, Vector<double> R, Vector<int> S, Vector<double> V, double t) {
			
			double dr = pw.s*0.977;
			
			for (int i=0; i<S.x; i++)
				for (int j=0; j<S.y; j++) 
					for (int k=0; k<S.z; k++){
					
						P pi;	
					
						pi.R.x = dr*i + (j%2==0 ? 0: 0.5*dr) + (k%2==0 ? 0: 0.5*dr);
						pi.R.y = 0.5*sqrt(3)*dr*j + (k%2==0 ? 0: sqrt(3.0)/6.0*dr);
						pi.R.z = sqrt(6.0)/3.0*dr*k;
					
						pi.V = Vector<double>(t*(1-2*rnd()), t*(1-2*rnd()), t*(1-2*rnd()));
					
						/*
						pi.R.x-= 0.5*S.x*dr;
						pi.R.y-= 0.5*S.y*0.5*sqrt(3)*dr;
						pi.R.z-= 0.5*S.z*sqrt(6.0)/3.0*dr;
						*/
						pi.R  += R;
						
						pi.V  += V;

						p.push_back(pi);	
					
				}
		}
		
		
		
		MyVisualizator(){
			eye.R.z = eye.RM.z = 30;
			
			
			pt.resize(1);
			
			pt[0].color = Color(1,0.2,0.1,1);
			pt[0].m = 1;
			
			pw[0][0].s  = 1;
			pw[0][0].e  = 1;
			pw[0][0].rc = 3;
			
			//packToHCP(p, pw[0][0], Vector<double>(1,1,1), Vector<int>(160,160,160),  Vector<double>(+5,0,0), 0);
			packToHCP(p, pw[0][0], Vector<double>(14,14,14), Vector<int>(32,32,32),  Vector<double>(-1,0,0), 0);
			
			//packToHCP(p, pw[0][0], Vector<double>(14,65,14), Vector<int>(20,20,20),  Vector<double>(0,-10,0), 0);
			/**/
			
			perspectiveF = false;
			scale = -1;
			pointSize = 2.5;
			
			integrator.parent = this;
			integrator.start();
			coolingF = false;
			
			
			grid_step = 1;
		};
		
		class Integrator : public QThread {
			public:
				
				MyVisualizator *parent;

				typedef long unsigned int luint;
				
				dVector w;
				bool exitFlag;
				bool ssF, eF, pF, firstTimeF;
//				double dt;
				double pot;
				double lU, eU, mU, tU;
				double etot, ek, ep, etot0, eerr;
				double temp;
				uint ntot;
				unsigned int step, nsteps;
				Timer tr, trf;
				unsigned int nstep;
				long unsigned int flops1, flops2, flops;
				
				Integrator():
				ssF(false),
				eF(false),
				pF(false),
				firstTimeF(true),
				nsteps(1<<32-1),
				flops1(0),flops2(0),flops(0)
				{
					lU = 1e-9; // m [1nm]
					eU = 1.60217646e-19; // J [1ev]
					mU = 1.660538921e-27; // kg [amu]
					tU = lU*sqrt(mU/eU); // sec
			
				}


			
			
			void calc_force(vector<P*> &pi) {
					
					
					
					for (int i=0;i<pi.size()-1;i++)
						for (int j=i+1;j<pi.size();j++) {
						
						
							double dx = pi[i]->R.x - pi[j]->R.x; // 1
							double dy = pi[i]->R.y - pi[j]->R.y; // 1
							double dz = pi[i]->R.z - pi[j]->R.z; // 1
								
							double r2 = dx*dx + dy*dy + dz*dz; // 5

							if (r2>drij2max) continue; 
							
							double r4 	= r2*r2; // 1
							double r6i 	= 1.0/(r2*r4); // 2
							double r12i = r6i*r6i;  // 1
								
							//pot += r12i - 2*r6i;  // 3
							double f = +e*r6i*r4*(r12i-r6i); // 4
							
							double dxf = dx*f;
							double dyf = dy*f;
							double dzf = dz*f;
							
							pi[i]->F.x += dxf; // 2
							pi[i]->F.y += dyf; // 2
							pi[i]->F.z += dzf; // 2
						
							pi[j]->F.x -= dxf; // 2
							pi[j]->F.y -= dyf; // 2
							pi[j]->F.z -= dzf; // 2
							
						}
				}


				void calc_force(vector<P*> &pi,  vector<P*> &pj) {
					
					//drij2max = SQR(3.0); 
					
					for (int i=0;i<pi.size();i++)
						for (int j=0;j<pj.size();j++) {
						
						
							double dx = pi[i]->R.x - pj[j]->R.x; // 1
							double dy = pi[i]->R.y - pj[j]->R.y; // 1
							double dz = pi[i]->R.z - pj[j]->R.z; // 1
								
							double r2 = dx*dx + dy*dy + dz*dz; // 5

							if (r2>drij2max) continue; 
							
							double r4 	= r2*r2; // 1
							double r6i 	= 1.0/(r2*r4); // 2
							double r12i = r6i*r6i;  // 1
								
							//pot += r12i - 2*r6i;  // 3
							double f = +e*r6i*r4*(r12i-r6i); // 4

								
							pi[i]->F.x += dx*f; // 2
							pi[i]->F.y += dy*f; // 2
							pi[i]->F.z += dz*f; // 2
						
							
						}
					
				}


				void run(){
					

					vector<PT> &pt = parent->pt;
					vector<P> &p = parent->p;
					
					for(int i=0; i<p.size(); i++) p[i].F*=0;

					exitFlag = false;
					
					step=0;
					
					ntot=p.size();
					nsteps = 15.0/dt;

					do {

						tr.set();
						
						
						for(int i=0; i<p.size(); i++) {
							p[i].R+= (p[i].V + .5*p[i].F*dt)*dt;
							p[i].V+= .5*p[i].F*dt;
						}
						
						// boundaries
						for(int i=0; i<p.size(); i++) {
							
							if (p[i].R.x<0) { p[i].R.x=0;  p[i].V.x*=-1; }
							if (p[i].R.y<0) { p[i].R.y=0;  p[i].V.y*=-1; }
							if (p[i].R.z<0) { p[i].R.z=0;  p[i].V.z*=-1; }
							
							if (gsizex*gstepx<p[i].R.x) { p[i].R.x=gsizex*gstepx;  p[i].V.x*=-1; }
							if (gsizey*gstepy<p[i].R.y) { p[i].R.y=gsizey*gstepy;  p[i].V.y*=-1; }
							if (gsizez*gstepz<p[i].R.z) { p[i].R.z=gsizez*gstepz;  p[i].V.z*=-1; }
							
						}
						
						
						for (int xi=0; xi<gsizex; xi++)
						for (int yi=0; yi<gsizey; yi++)
						for (int zi=0; zi<gsizez; zi++)
							g[xi][yi][zi].resize(0);
						
						for(int i=0; i<p.size(); i++)
							g[uint(p[i].R.x/gstepx)][uint(p[i].R.y/gstepy)][uint(p[i].R.z/gstepz)].push_back(&p[i]);
						
						trf.set();
						
						for (int xi=0; xi<gsizex; xi++)
						for (int yi=0; yi<gsizey; yi++)
						for (int zi=0; zi<gsizez; zi++) {
							
							
							vector<P*> &pi = g[xi][yi][zi];
							
							if (pi.size()==0) continue;
								
							for (int i=0; i<pi.size(); i++) pi[i]->F*=0;
							
							calc_force(pi);

							
							for (int xj= (xi==0 ? xi:xi-1); xj<=(xi<gsizex-1 ? xi+1: gsizex-1); xj++)
							for (int yj= (yi==0 ? yi:yi-1); yj<=(yi<gsizey-1 ? yi+1: gsizey-1); yj++)
							for (int zj= (zi==0 ? zi:zi-1); zj<=(zi<gsizez-1 ? zi+1: gsizez-1); zj++) 
							if (xi!=xj || yi!=yj || zi!=zj )
							calc_force(pi,  g[xj][yj][zj]);
							
						}
						trf.get();
						
						for(int i=0; i<p.size(); i++)
							p[i].V+= .5*p[i].F*dt;
						
						
						step++;

						/*
						if (step%uint(15.0/dt)==0) {
							for (int i=0; i<p.size(); i++)
								printf("%.4e %.4e %.4e\n", p[i].R.x, p[i].R.y, p[i].R.z);
						
						}
						*/	
						tr.get();
						
						//printf("%9i %10.4f %10.4f\n", ntot, trf.sin, tr.sin);
						
					} while (step<nsteps && !exitFlag);
					parent->close();
				}
				
				
				~Integrator(){ exitFlag = true; }
				
		} integrator;
		
		
		virtual void printInfo(){
			
			Visualizator::printInfo();
			
			
			char str[512];
			sprintf(str, 
			"step: % .2e\ndt: %.2e\nt: %.2f\nsteps/s: %g\nsteps/s: %g (force)\n\np.size(): %i", 
					  (double)integrator.step, dt, integrator.step*dt,
					  1.0/integrator.tr.sin, 1.0/integrator.trf.sin, integrator.ntot		  
					  );
			
			glColor3f (.4,.4,.4);
			glPrint(10,85, str, GLUT_BITMAP_HELVETICA_12);
			
		}
		
		void paintGL() {

			glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
			
			Visualizator::paintGL();
				
			uint k=0;
			glBindBuffer(GL_ARRAY_BUFFER_ARB, k+1);
			glBufferData(GL_ARRAY_BUFFER_ARB, p.size()*sizeof(p[0]), &(p[0].R), GL_STATIC_DRAW_ARB);
			drawPointArray(k+1, p.size(), pt[k].color, pt[k].color, pointSize, GL_ONE, GL_ONE, false, program, texture, sizeof(p[0]), 0, GL_DOUBLE);

		};
		
		void keyPressEvent(QKeyEvent *event){
			Visualizator::keyPressEvent(event);
			
			if (event->key() == Qt::Key_PageUp)	scale+= -1.0/scale_steps;
			if (event->key() == Qt::Key_PageDown) scale+= +1.0/scale_steps;
			
			if (event->key() == Qt::Key_C) coolingF=!coolingF;
			
			
		}
		
		
		virtual void drawGrid(){
			
			glColor4fv(grid_color);
			
			
			for (int i=0; i<=gsizex; i++) {
				
				glBegin(GL_LINE_LOOP);	
				
				glVertex3f(gstepx*i,0, 0);
				glVertex3f(gstepx*i,gstepy*gsizey,0);
				glVertex3f(gstepx*i,gstepy*gsizey,gstepz*gsizez);
				glVertex3f(gstepx*i,0,gstepz*gsizez);
				
				glEnd();	
			}
			
			
			for (int i=0; i<=gsizey; i++) {
				
				glBegin(GL_LINE_LOOP);	
				
				glVertex3f(0, gstepy*i, 0);
				glVertex3f(gstepx*gsizex, gstepy*i,0);
				glVertex3f(gstepx*gsizex, gstepy*i,gstepz*gsizez);
				glVertex3f(0,gstepy*i,gstepz*gsizez);
				
				glEnd();	
			}
			
			
			for (int i=0; i<=gsizez; i++) {
				
				glBegin(GL_LINE_LOOP);	
				
				glVertex3f(0, 0, gstepz*i);
				glVertex3f(gstepx*gsizex,0,  gstepz*i);
				glVertex3f(gstepx*gsizex, gstepy*gsizey, gstepz*i);
				glVertex3f(0,gstepy*gsizey, gstepz*i);
				
				glEnd();	
			}

	
			
		}
	
};




int main(int argc, char *argv[]) {
	
	
	
	
	
	
	
	glutInit(&argc, argv); //initialize glut
	QApplication app(argc, argv);
	
	MyVisualizator visualizator;
	
	bool windowF(true);
	visualizator.argv=argv;
	
	for (int i=1;i<argc;i++)
		if (!strcmp(argv[i],"-nw")) windowF=false;
		else if (!strcmp(argv[i],"-ssF")) { visualizator.integrator.ssF=true;  mkdir("ss", 0777);  }
		else if (!strcmp(argv[i],"-nsteps")) visualizator.integrator.nsteps = atoi(argv[++i]);
		else if (!strcmp(argv[i],"-gstepx")) gstepx = atof(argv[++i]);
	
	visualizator.resize(800, 800);
	visualizator.setWindowTitle( QString("Visualizator Lite") );
	if (windowF) visualizator.show();
	
	return app.exec();
	
}
